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ABSTRACT 

The extrasolar planets discovered to date possess unexpected orbital elements. Most 
orbit their host stars with larger eccentricities and smaller semi-major axes than simi- 
larly sized planets in our own solar system do. It is generally agreed that the interaction 
between giant planets and circumstellar disks (Type II migration) drives these planets 
inward to small radii, but the effect of these same disks on orbital eccentricity, e, is 
controversial. Several recent analytic calculations suggest that disk-planet interactions 
can excite eccentricity, while numerical studies generally produce eccentricity damping. 
This paper addresses this controversy using a quasi-analytic approach, drawing on sev- 
eral preceding analytic studies. This work refines the current treatment of eccentricity 
evolution by removing several approximations from the calculation of disk torques. We 
encounter neither uniform damping nor uniform excitation of orbital eccentricity, but 
rather a function de/dt that varies in both sign and magnitude depending on eccentricity 
and other solar system properties. Most significantly, we find that for every combination 
of disk and planet properties investigated herein, corotation torques produce negative 
values of de/dt for some range in e within the interval [0.1, 0.5]. If corotation torques 
are saturated, this region of eccentricity damping disappears, and excitation occurs on 
a short timescale of less than 0.08 Myr. Thus, our study does not produce eccentricity 
excitation on a timescale of a few Myr - we obtain either eccentricity excitation on a 
short time scale, or eccentricity damping on a longer time scale. Finally, we discuss 
the implications of this result for producing the observed range in extrasolar planet 
eccentricity. 
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1. Introduction 



Our solar system contained the only known planets orbiting a main sequence star until 1995, 
when a planet was found orbiting the Sun-like star 51 Peg via a periodic Doppler shift in the 
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star's spectrum. The discovery of more than 150 additional extrasolar planets has overturned our 
understanding of what constitutes a typical planetary system. The extrasolar planets discovered to 
date are Jupiter-sized, possess a wide range of orbital eccentricities (0 < e < 0.9), and orbit their 
host stars with small semi-major axes (0.03 AU < a < 6 AU). In contrast, similarly sized planets 
in our solar system (Jupiter and Saturn) live in nearly circular orbits at 5 and 10 AU. 

These discoveries prompted corresponding shifts in solar system evolution theory. We previ- 
ously believed that planets formed in, or near, their current orbits. However, ice, which plays a 
significant role in forming the dense solid core of a gas giant, will not condense within 3 AU of 
a Sun- like star, which in turn implies that gas giant planets are unlikely to form with the small 
semi-major axes they possess. We now believe that massive planets form outside of this "snow 
line" and subsequently move inward via interactions with a circumstellar disk, a process known as 
planetary migration. There are two limiting cases in migration theory in a laminar disk: Type I 
migration, in which a planet lacks sufficient mass to clear a gap in the disk material and is driven 
inward by a density wake in the disk, and Type II migration, in which massive planets do clear 
a gap and are driven inward by resonances between the planet and material in the remainder of 
the disk. If the disk is not laminar. Type I migration can be overwhelmed by turbulent effects 
(Laughlin et al. 2004, Nelson 2005). 

In the standard description of planetary migration, rings of disk material in resonance with a 
planet in a nearly circular orbit exert torques on the planet, driving it inward and further decreasing 
its eccentricity e (Goldreich and Tremaine 1980, hereafter GT80). This scenario, while appropriate 
for our solar system, never allows the planet to possess large orbital eccentricity, and hence does not 
explain the large eccentricities of many detected planetary orbits. As a result, one of the current 
challenges in planetary migration theory is to produce the wide range in orbital eccentricity observed 
in the extrasolar planet population. The observed range in eccentricity can be explained by disk 
torques acting in conjunction with interactions with a second planet (Moorhead and Adams 2005, 
hereafter MA05; see also Rasio and Ford 1996, Thommes and Lissauer 2003, Adams and Laughlin 
2003). Here we investigate whether eccentricity excitation can take place when we remove the small 
eccentricity assumption from the GT80 formalism. 

The combination of planet-planet scattering and simplified disk effects (constant damping 
rates of semi-major axis and eccentricity) successfully reproduces the observed a-e distribution of 
the extrasolar planets (MA05). The disk drives both planets inward while interactions between 
the planets compete with the disk's eccentricity damping; the result is both small a and a large 
range in e. In treating the disk as a source of constant eccentricity damping, this past work was 
based on a host of analytic (e.g., GT80) and numerical (Trilling et al. 1998, Nelson et al. 2000, 
Papaloizou et al. 2001) studies. However, several recent studies (Goldreich and Sari 2003, hereafter 
GS03; Ogilvie and Lubow 2003, hereafter OL03) suggest that disk torques may excite rather than 
damp eccentricity. In this paper, we investigate under what conditions nearly Keplerian disks - 
like those observed in nearby regions of star formation - will give rise to eccentricity excitation or 
damping. We then outline a method for incorporating the results into numerical simulations in 
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order to investigate to what extent the a-e distribution changes under this modification. 

This paper combines the formahsm of three key analytic papers on planet migration: GT80, 
GS03, and OL03. Specifically, we obtain our torque formulas from GT80, the form of the eccentricity 
derivative from GS03, and our treatment of saturated corotation resonances from OL03. We extend 
this combined treatment primarily by performing a full calculation of the coefficients of the cosine 
expansion of the planet's perturbation to the overall gravitational potential, (j^^jy^iP), assuming 
neither small eccentricity nor large azimuthal wavenumber (as in past studies). We find that the 
function has a complicated dependence on eccentricity (see section 2). Each resonant torque 

depends on 4>^^{I3)., and as a result, the eccentricity time derivative is also a complicated function of 
eccentricity; in particular, de/dt attains both positive and negative values depending on the current 
value of eccentricity. Several properties of the planet's cleared gap, including its width, placement 
around the planet, and the degree to which it is cleared, alter the shape of de/dt. In this paper, 
we present eccentricity time derivatives for a variety of gap widths and shapes, seeking common 
behaviors. Additionally, if corotation resonances become saturated, we find de/dt to be almost 
exclusively positive. In this manner, we attain a greater understanding of how the production 
of either eccentricity excitation or damping depends on current planet and disk properties. The 
paper concludes with a discussion of this finding in the context of recent analytic studies as well as 
numerical studies, which generally produce eccentricity damping. 



2. Methods and Initial Conditions 

2.1. Disk Properties 

For the sake of definiteness, we consider a IMj planet orbiting a 1 Mq star, embedded in 
a 0.05 Mq circumstellar disk with radius 30 AU. We assume that the unperturbed disk surface 
profile density falls off with distance from the central star as a power law, S = So(r/ro)-^/2. We 
choose this power-law index for comparison with the many numerical studies that do the same; 
additionally, the results are largely insensitive to this choice. We can easily verify that the small 
chosen disk mass does not disturb the overall Keplerian rotation curve ^{r), i.e., the condition 

W-r « ' 

is satisfied; the largest value GS/fi^r attains for any value of r in our given range is 0.04. 

We assume that the radial temperature profile also obeys a power law, T = Toir /ro)~^/^ , 
where To = 50 K at the snow-line, tq = 7 AU. Finally, we specify the properties of the disk 
material, the standard viscosity parameter for accretion disks a = 10"'^ and the ratio of specific 
heats 7 = 1.4. These parameters determine the extent to which planets are able to clear gaps in 
the disk and the strength of the torques produced by disk-planet resonances (see GT80, Shu 1992, 
Lin and Papaloizou 1993). 
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We assume a flat (cold) disk in our analysis. In addition, our torque formulas require that the 
disk satisfy the more stringent condition, 

h (cs/n) 
r r 

for each resonance of order m, where h is the disk scale height and Cg is the sound speed. The 
largest value h/r attains is 0.07, where h/r ~ 0.04(r/l AU)"^/^. A 1 Mj planet at 1 AU allows for 
only 5-10 resonances of each type, and so we can safely use the flat disk torque formulas. Notice, 
however, that smaller planets produce smaller gaps, which in turn require larger values of m. As a 
result, this treatment cannot be extrapolated to arbitrarily small planet masses. 



2.2. Disk Torques 

For a Keplerian disk, the eccentricity evolution due to disk torques is given by (GS03) 

dt Mpy/GM^a ^' ^ ' 

where a, e, and Mp are the semi-major axis, eccentricity, and mass of the planet, m is the azimuthal 
wavenumber of the pattern speed (see below), and To is the portion of the torque exerted by the disk 
on the planet that corresponds to the Fourier component of the planet's potential with azimuthal 
wavenumber m and pattern speed {l/m)Qp. 

The system contains an infinite number of resonances, and corresponding torques, between the 
disk and the planet. For example, if we assign these resonances wavenumbers (m, i), the location of 
a corotation resonance is given byre = a[m/(£)]3/2. r^y^^ 

massive planets of interest here will clear 
large gaps in the disk, providing, for each value of £, an upper limit on the number of resonances 
contributing to the total torque. The shape of the gap also affects the eccentricity evolution; 
Lindblad resonances are proportional to surface density, and corotation resonances are proportional 
to the radial derivative of surface density. We first calculate the location and torque of each 
resonance (and discuss the perturbing function's coefficients <p^^{P) in detail). We then present 
our method for translating a given gap profile into an upper limit on the number of contributing 
resonances. 

Resonances occur, and torques are exerted, where the motion of a ring in the disk matches the 
pattern speed ^£,m of the planet. For a Keplerian disk, this condition takes the form 

= i}p + (i — m)Kp/m = — Qp . (4) 

m 

We must consider both Lindblad and corotation resonances in our calculations, as described above. 
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2.3. Lindblad Resonances 



Lindblad resonances occur at radii r where it K{r)/m 
Keplerian disk, this condition can be written in the form 



£,nii 



nir) 



m 



mil 



n 



£,m ■ 



where m > 0. For a 
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The radius of a Lindblad resonance is denoted ri. In the above equation and throughout this 
discussion, we take the top sign for an outer Lindblad resonance and the lower sign for an inner 
Lindblad resonance. The expression for the disk torque due to a Lindblad resonance in a cold, 
Keplerian, non-gravitating disk is given by 

\ 2" 



rmr 



3(1 ±m) 



dr 



£,m, „ .p 

=F 2m(pp 



(6) 



rL 



where (p^^ is the {I, m) component of the cosine expansion of the disturbing potential produced 
by the planet (see GT80), and is discussed in further detail below. 



2.4. Corotation Resonances 



Corotation resonances occur at radii r where 



(7) 



We denote the radius at which we encounter a corotation resonance rc- The expression for the disk 
torque due to a corotation resonance (GT80) in a cold, Keplerian, non-gravitating disk is given by 



4m7r 



d /S 



n{r) dr \n 



(8) 



rc 



2.5. Determining the Expansion Coefficients 

Just as the disk is not massive enough to be self-gravitating, the orbiting planet is not massive 
enough to alter the 1/r potential produced by the central star. As a result, we treat the planet's 
gravitational potential as a small perturbation in the overall potential, and consider the cosine 
expansion of the potential. The elements of this expansion correspond to resonances between the 
planet and disk material at different radii. 

The disturbing potential cjF produced by an orbiting planet moving in the plane of the disk is 
well known (for example, Murray and Dermott 2001), and is given by 

(^P(r,0,t) = GMp(^-^-^^-\ (9) 
\ 7. J |r _ rp| y 
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GMp 



rp cos {9 — 6p) 



/r'^ + rp — 2rrp cos {9 — Bp) 
We can expand in terms corresponding to pattern speeds 0^^: 

oo oo 

</.P(r,e,0= Yl ^Uir) cos {m9-mpt). (10) 



£=—oo m=0 

Defining /3 = r/a, we can write the expansion coefficients in the form 



GMp 



2n 1-2-K 



COS (m9 — £{S^ — esin^))(l — ecos.^) 
y^/3^ + {l-ecosO^-2Pg{e,0/ 



(11) 



5(^)0 = (cos ^ — e) cos (0) + (Vl-e2 sin^) sin(e) 



Note that P is not a free variable, but is determined by m, i, and the type of resonance. As 
a result, 4>im(f^) ^ non-trivial function of m, i, e and the type of resonance, as both the planet 
mass Mp and the semi-major axis a are prefactors. 

Calculating the coefficients (p^^ requires the evaluation of an oscillatory two-dimensional in- 
tegral. One method of sidestepping this problem is to expand <^^„(r) to first order in e (GT80). 
However, extrasolar planets often display large orbital eccentricity, and one goal of migration the- 
ory is to explain these large eccentricities. It is therefore important to consider the validity of the 
linear approximation of 4>£j^ for the full range of eccentricity e. Figure 1 displays the behavior of 
the numerically determined and approximate expressions for one particular coefficient, i;^4 3(riLR), 
as a function of eccentricity. Clearly, an approximation to first order in eccentricity is only useful 
for eccentricities e < 0.3. In the context of GT80, this approximation was acceptable as researchers 
were mainly interested in explaining the properties of our own solar system, where planets have 
nearly circular orbits. However, many extrasolar planets possess large eccentricities, and thus the 
small e approximation necessarily breaks down. 

An alternative method for avoiding the full calculation of is to take the large m limit, 
where the Lindblad and corotation resonances have the same functional dependence on (t>£^{r) 
(see Eqs. 6 and 8). In this limit, one can evaluate the degree to which Lindblad and corotation 
resonances cancel or add (GS03). The torque is proportional to m, so large m resonances are more 
significant. This approximation works well for less massive planets, which clear small gaps and 
allow large m resonances to contribute to the torque. On the other hand, a 10 Mj planet clears a 
wide enough gap that only a few resonances exist, and even a 1 Mj planet only allows for only 5 
to 10 resonances of each type. Thus, an important part of this study is to calculate the coefficients 
0^^(r) to much higher order accuracy for larger values of eccentricity, e. 
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To conclude, for the currently observed sample of extrasolar planets, making low e or large 
m approximations may not produce the correct eccentricity evolution for any given system; on 
the other hand, calculating the exact solution for (j}£^{r) is often not computationally feasible. 
Therefore, we embark on a semi-analytic investigation of the eccentricity evolution for intermediate 
values of eccentricity and for a number of representative disk-planet configurations. We look for an 
overall pattern, both to understand eccentricity evolution and to incorporate into future calculations 
of planetary migration. 



2.6. The shape of the function 4>^i^{li) and its radial derivative 

Here we note several characteristics of and its derivative as a function of e, as shown 

in Fig. 1. The overall shape consists of three distinct regions: a nearly linear decrease with e for 
small e, a cusp located at an intermediate value e = eo, and a smooth oscillation for larger values 
of e. As m increases, the location of the cusp moves to the left, and the number of oscillations 
increases roughly in proportion with m. The cusp at eo = 1 — /? is of particular interest because 
the derivative dcp^ / d(3 becomes unbounded at e = eo. 

The function <^£^(/3) is bounded, and we find, as expected (GT80) that the overall amplitude 
of the resonance falls off as \i — ml increases (GT80). However, the derivative dcf)^ ^{(3) / dr is 
unbounded at the location of the cusp in 4>^^{f}), eo. Furthermore, = \1 — f3\ = \1 — {m — 
(for inner Lindblad resonances, for example), and so is dense like the rational numbers. In other 
words, 

Ve, /iGM,{</.P : \dAl^{e) 



>h}^0. (12) 



Using this information alone, it is incorrect to ignore higher values of |/ — m|. However, any 
physical disk will exhibit a smooth response to the planet's stimulus. To mimic this effect, we 
manually smooth the potential (and therefore the cusp in </<£„) over the largest relevant physical 
scale. The viscous length scale is given by 

^"^[mi-dn/dr]) ' ^^^^ 

where v is the characteristic kinematic viscosity (OL03). For m = 1 and r = 1 AU, (5jy = 0.004, 
and for any m, r of interest, 5i, < 0.01. On the other hand, the size of the disk scale height is 
h/r ~ 0.04. We will present results using first the viscous length scale as our chosen smoothing 
length, and then using the disk scale height. This is to demonstrate the effect that increasing the 
smoothing length has on our results. Smoothing over these physical scales prevents the derivative 
from attaining arbitrarily large values, and higher order resonances may be safely discarded. 



Furthermore, we may apply a low order correction for our treatment of the disk as two dimen- 
sional by averaging over the disk scale height (Menou and Goodman, 2004). The disk scale height 



- 8 - 



is about four times as large as the viscous scale length, and therefore we will discard the viscous 
scale length and smooth all of our functions over the disk scale height. 

It is important to note that the amplitude of the resonances decreases much more rapidly for 
small e than for large e. Thus, it is important to include resonances for which |^ — m| > 1 in order to 
calculate de/dt accurately. However, the shape of 4>^m{P) does not change rapidly with |£ — m| (see 
Fig. 2), and higher order resonances tend to add to a similar shape as do lower order resonances. 
Therefore, we take the approach of using |£ — m| < 1 resonances to determine the shape of de/dt, 
and then calculating de/dt more accurately for key values of e using all contributing i, m. 

3. Results 

Using the framework developed in the above section, we can calculate the eccentricity time 
derivative for a given Jovian planet in a given thin disk. The most important parameters in this 
calculation are the mass, semi-major axis, and eccentricity of the planet, and the surface density 
profile in the vicinity of the planet's cleared gap. Our algorithm can be summarized as follows: [1] 
We define the radial surface density profile of the disk and the orbital properties of the planet. [2] 
We evaluate the formulae for (j)^^ and its derivative (using the Gauss-Kronrod numerical integration 
function in Mathematica), the torques, and the eccentricity derivative resulting from each resonance, 
smoothing over the disk scale height to partially account for the disk's finite thickness (Menou 
and Goodman 2004). [3] We sum the contributions to the total eccentricity derivative over all 
resonance types and values of m. For each set of input parameters (disk properties plus orbital 
parameters of the planet), we obtain de/dt as a function of e. Issue may be taken with any individual 
surface density profile we present in this section. However, we present de/dt for a variety of gap 
architectures, noting shared characteristics. 

3.1. A Sharp-Edged Gap 

We first present the simplest case: a 1 Mj planet orbiting at a = 1 AU in a gap with sharply 
defined edges. The strongest resonances occur for values of i = m or i = mzizl, limiting the number 
of values of i we must consider (GT80). The azimuthal wavenumber m may still take on any value, 
but, as m increases, the location of the resonance approaches the semi-major axis of the perturbing 
planet. When the planet clears a clean gap, the edges of the gap place a physical upper limit on 
the number of values of m that must be considered. The resulting function de/dt versus e (shown 
in Figs. 3 and 4) is the result of summing the contributions of five to ten resonances of each type 
(where type refers to, for example, an ^ = jtt, -|- 1 inner Lindblad resonance). 

In Fig. 3, de/dt is not smoothed over the disk scale height, but is instead smoothed over 
the viscous length scale 5^ ~ 0.01. Due to the small smoothing length, the features of the (j)^^ 
functions are clearly visible in Fig. 3; the peaks correspond to the singularities in dcj)^ ^{(3) / d(3 
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smoothed over the viscous length scale, and the smooth right half of the plot is produced by the 
combination of the smoothly oscillating portions of the and d4>^ ^{(3) / df3 functions. In Fig. 

4, the features of the functions are smoothed out over the disk scale height, h/r ~ 0.04, and 
thus not clearly visible in the resulting function de/dt. However, the peaks and troughs in de/dt for 
e < 0.5 generally correspond to the singularities in dcf)^ ^{fi) / dp , and peaks and troughs in de/dt for 
e > 0.5 correspond to the smoothly oscillating portions of the </'^„ (/3) and dcl)^ ^ {P) / d(3 functions. 
The most important result obtained from this plot is that planet-disk interactions do not result 
in uniform eccentricity damping or excitation. The behavior of de/dt is complicated, leading to 
varying levels of damping or excitation depending on the current value of e. 

There exist sizable bodies of both analytic and numerical work on the subject of planet-induced 
gaps in circumstellar disks. In the context of this study, these works translate into a wide array 
of possible inputs for the surface density in the vicinity of the planet. A planet on an eccentric 
orbit is likely to have a wider gap than a planet in a circular orbit, as is assumed in many of these 
studies, and so we also perform a calculation of de/dt for an artificially widened gap. The shape of 
the surface density profile has two important effects on the resulting eccentricity evolution equation 
de/dt: [1] Since each torque depends on either the surface density or its radial derivative, altering 
the shape of the surface density profile will (linearly) change the relative heights of the features in 
the function de/dt. At the same time, if the width of the gap is unchanged, the overall pattern 
will also be largely unchanged. [2] If we instead alter the gap width, the shape of the function 
de/dt changes as resonances are included or excluded. To understand the general implications of 
gap width for eccentricity evolution, we present gap architectures with varying gap widths. 

The result of slightly narrowing the cleared gap - a change in the sign of de/dt - is also 
displayed in Fig. 4. However, we are most interested in calculating de/dt for intermediate values 
of eccentricity (and this is the region where our assumptions are most accurate), so this change is 
not significant within the context of our calculations. 

3.2. The EfTect of Residue in the Gap 

When some material remains within the gap, the radius at which disk material no longer exerts 
torques on the planet is unclear. To investigate eccentricity evolution in such a gap, we choose a 
gap shape resembling that produced by the numerical calculations of Bate et al. (2003), where 
residual material remains within the gap (see Fig. 5). We then raise the limit on the number of 
contributing resonances (by raising the allowed values of m) until de/dt converges to a single value. 

As mentioned, the radius at which disk material no longer exerts torques on the planet is un- 
clear. Therefore, we present de/dt with and without torques resulting from coorbital disk material. 
Figure 6 displays [de/dt)/e as a function of e where only non-coorbital resonances are allowed to 
contribute. We find that {de/dt)/e does not converge as e approaches zero, but instead approaches 
+00. When we include the effects of coorbital resonances (Fig. 7), de/dt switches sign for small 
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e, resulting in eccentricity damping rather than excitation. Notice, however, that in both cases, 
eccentricity damping results for all values of eccentricity in the range 0.1 < e < 0.8. 

3.3. Other configurations 

As mentioned previously, the overall shape of the function de/dt versus e is fairly uniform 
for all surface density profiles that contain a sharp-edged gap of a given width centered on the 
planet. However, there are a number of scenarios in which a planet may lie off-center in the gap: 
accretion may clear the inner disk, photoevaporation or a close encounter with another star may 
strip the outer disk, or multiple planets may produce overlapping gaps in the disk. Of these three 
possibilities, the situation in which two planets have overlapping gaps is of the greatest interest 
because two planets orbiting in a close resonance are certain to produce a gap in which neither 
planet lies at the center. In contrast, it requires more fine-tuning for other processes to deplete the 
disk from one of its edges exactly to the semi-major axis of the planet. 

To study the effect of a two-planet gap, we place two 1 Mj planets in orbits at ai =1 AU 
and a2 =1.59 AU, and thus in a 2-1 mean motion resonance. The resulting gap in disk material is 
displayed in Fig. 8, and the eccentricity evolution it produces is displayed in Fig. 9. We find that 
the gap asymmetry does not affect the planets equally; while the inner planet experiences a drastic 
change in both the shape and magnitude of its de/dt function, the outer planet's eccentricity evolves 
in much the same way a lone planet's eccentricity evolves. We find, upon closer inspection, that 
this result is due to the fact that the contribution to de/dt from the m = 2, ^ = 1 corotation and 
Lindblad resonances are significantly greater than any other resonance contributions for a single 
planet, and that these resonances are located on the outer side of the gap. For the inner planet in 
a 2-1 resonance, there is no material at this location, but for the outer planet, the location of the 
2-1 resonance lies within the disk material. We can deduce from this that, for any planet orbiting 
in a disk, a particularly strong resonance exists between the planet and a ring of disk material 
inside the planet's orbit. It is then important to combine this effect with that of the mean motion 
resonance (Murray and Dermott 2001, Thommes and Lissauer 2003, Lee and Peale 2002) between 
the two planets when performing semi-analytic studies of solar system evolution such as that in 
MA05. 

If these results are to be used in future numerical calculations such as MA05 it is important to 
understand how the eccentricity evolution depends on orbital parameters. By completing a large 
range of calculations, we have determined the following: [1] As the semi-major axis a decreases, 
the surface density in the vicinity of the planet increases, and the planet sits in a deeper, narrower 
trough. The function de/dt versus e obtains an increased number of peaks and troughs, smoothed 
over the disk scale height, in the small e region as we include more resonances, and an increased 
magnitude of de/dt everywhere due to the larger surface density at small radii. As a increases 
rather than decreases, these effects are reversed. [2] As the planet mass decreases, the planet clears 
a narrower and narrower gap. Once again, the function de/dt versus e obtains an increased number 



- 11 - 



of peaks and troughs, again smoothed over the disk scale height, and experiences an decreased 
overall magnitude due to the dependence de/dt oc Mp. (This dependence arises from Eqs. 3 and 
11; In Eq. 3, the formula de/dt explicitly contains the term 1/Mp, and we see from Eq. 11 that 
[4>\rnf' contributes M|,.) 

3.4. Saturation of Corotation Resonances 

The torque of the planet on resonances in the disk deposits or removes angular momentum from 
the resonance locations in the disk. Lindblad resonances lose this change in angular momentum 
through a wave flux, but corotation resonances do not. As a result, they can become saturated and 
thereby exert a reduced torque on the planet. In this section, we assume that coorbital corotation 
resonances will be completely saturated, and thus provide no contribution to the overall eccentricity 
evolution (Balmforth and Korycansky, 2001). We account for the effects of corotation saturation 
using the prescription of OL03. 

The degree of saturation for a particular (£, m) non-coorbital corotation resonance can be 
expressed in terms of a single parameter p (OL03), which is given by 

_ ( -d\nr \ 4>lm{li) ( [-dn/dr] \ 

V ) k2 \ u/m ) ' ^ ' 

where the viscosity v is related to the parameter a by the standard Shakura-Sunyaev formulation, 
V = aO.h?' . Therefore, assuming the disk is Keplerian, and that the disk scale height h = Cg/^, the 
expression for p becomes 

2 p 

Using the quantities specified previously, we obtain p k, 0.284 4>^^{I3). Our sample calculation of 
04 3(^ilr) has an amplitude of about 5 (</>4 3(rc) has a similar amplitude), and so the saturation 
parameter p attains values of order unity for some values of e. Thus, it is important to investigate 
the effect that saturation may have on our results. 

The torque exerted in the corotation region is reduced by a factor /(p), where f{p) ~ 0.4019p~^/^ 
for p » 1, and /(p) w 1 — 2.044^^ for p ^1 (OL03). We interpolate between these two limits in a 
manner similar to that in GS03 to obtain 

fij)) = (1 + 0.3529/)5/V(l + 1.022p2)2 . (16) 

The result of including saturation effects is shown in Figs. 10 and 11. Thus, the calculated satura- 
tion level of corotation resonances results in eccentricity excitation over the full range of eccentricity. 
If coorbital resonances are allowed to contribute, a small interval of eccentricity damping exists for 
e < 0.05. On a side note, we learn that the substantial eccentricity damping in the low e limit in 
Fig. 11 is due to coorbital Lindblad resonances and not coorbital corotation resonances. 



2>mH 



ajkBTo V GM-. 



(15) 
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Taken at face value, our results imply that all single planets starting with e > 0.1 will be 
ejected from their systems or accreted on time scales shorter than 0.15 Myr. Furthermore, uniform 
eccentricity excitation would lead to no multiple planet systems (MA05), in contrast with the 
observed population of extrasolar planets which contains about 10-20% multiple systems. As a 
result, it is unlikely that corotation resonances are saturated to the degree described above. The 
determination of the degree of saturation is thus an important issue for the future. 

3.5. Eccentricity Damping at e = 0.3 

As mentioned earlier in this work, our goal is to compute the eccentricity time derivative for a 
variety of orbital parameters and gap shapes, looking for a common pattern. We find that for each 
set of inputs, the eccentricity time derivative for a single planet without saturation of corotation 
resonances is negative for some range in eccentricity between 0.1 and 0.5, and is always negative 
for e = 0.3. On the other hand, the eccentricity derivative for the same set of parameters, including 
e = 0.3, is always large and positive when corotation resonances become saturated. 

At this point in our study, we narrow our focus to e = 0.3 only. We now include resonances in 
our calculation for which \£ — m\ > 1 and recalculate de/dt for our different gap shapes. The result 
is displayed in Table 1. The corresponding gaps in surface density are displayed in Fig. 12. We see 
that de/dt maintains its sign while adopting a greater magnitude in all cases, as expected. 

It is clear that, using our approach and assumptions, we obtain, for intermediate eccentricity, 
either eccentricity damping or extremely rapid eccentricity excitation. Eccentricity damping at 
e = 0.3 will prevent planets from attaining the full range of eccentricity values that the observed 
sample of extrasolar planets possesses. On the other hand, if corotation saturation occurs according 
to the formula of Ogilvie and Lubow (2003), planets would attain very large eccentricities in a few 
thousand years, and thus would likely be ejected from their host systems within the Myr lifetime of 
the disk. Either some additional component (such as vorticity of the gap edges) must be included in 
calculations of eccentricity evolution, or some additional physical phenomenon, such as interactions 
between planets in multiple planet systems, must take place. 

4. Summary and Discussion 

In this paper, we have re-examined the manner in which circumstellar disks exert torques on 
giant planets and thereby change their orbital eccentricity. We have shown that both the small e and 
large m approximations for calculating eccentricity evolution are invalid for studying the current 
population of observed extrasolar planets. We have removed both of these approximations from 
the formalism and completed a full calculation of <p£m{P) (see Fig. 1). The resulting function de/dt 
attains both positive and negative values (see Figs. 3 and 4), reflecting its nature as a composition 
of </'£m(/3) functions. Thus, Type II migration produces neither uniform eccentricity excitation 



- 13 - 



nor uniform eccentricity damping, but can produce either excitation or damping depending on the 
combination of disk and planet properties. We have found that de/dt is influenced by gap properties, 
in particular gap width. As the gap width changes, the number of included torques changes and 
we see corresponding addition or subtraction of spikes in the de/dt function. The placement of the 
planet within the gap is also important. A planet in a 2-1 resonance with a second planet and 
thus far offset from the gap center, for instance, may experience substantially different eccentricity 
evolution than a planet roughly at the center of the gap. 

A gap with sharp edges places an upper limit on the number of contributing resonances. The 
exact width of such a gap strongly affects the eccentricity time derivative in the low eccentricity 
limit. For instance, a small change in the gap width produces eccentricity excitation rather than 
damping for small eccentricity (Fig. 4). Therefore, it is important to understand the shape of 
the planet's cleared gap to predict the eccentricity evolution in the small eccentricity regime using 
our approach. If the gap contains residual contributing material (Figs. 6 and 7) or if the gap is 
generated by two, rather than one, planets (Fig. 9), the shape of de/dt versus e changes drastically. 

This work has several limitations. First, we have assumed throughout that the disk is infinites- 
imally thin. The accuracy of our method relies on the inclusion of relatively few disk torques. As m 
increases, the flat-disk approximation for disk torques is less accurate. This effect is of concern in 
our analyses when e < 0.1. Second, our calculations converge extremely slowly for large eccentric- 
ity. Furthermore, while D'Angelo et al (2006) showed that the primary effect of moderate planet 
eccentricity was to widen the gap in the disk, the effects of a very eccentric planet (e > 0.5) on 
disk geometry are uncertain. With these limitations in mind, we present again our de/dt plots with 
data shown only for 0.1 < e < 0.5 (Fig. 13). This plot demonstrates that, as long as corotation 
resonances are unsaturated, eccentricity damping necessarily occurs for 0.2 < e < 0.5. There- 
fore, disk-planet interactions alone cannot produce the full observed range in extrasolar planet 
eccentricities, assuming corotation resonances are unsaturated. 

If, on the other hand, we find that if corotation resonances are saturated, we obtain almost 
exclusive eccentricity excitation (Figs. 10 and 11). However, the resulting excitation is so severe 
that a planet could only remain in orbit about its host star for less than a tenth of a million years 
(we obtain this timescale when we include all values of i, m in our de/dt calculation). Assuming 
that disks last for of order one million years, this scenario fails to explain how planets survive 
disk-planet interactions. 

There are several possible refinements and extensions of this work. In our treatment of the 
disk as strictly Keplerian, we neglect the effects of vortices near the gap edges. A careful treatment 
of vorticity may alter the strength of the corotation resonances and lead different values of the 
eccentricity derivative. Furthermore, it may be informative to extend the range in eccentricity for 
which this approach is valid, as the observed range in eccentricity for the extrasolar planets ranges 
from to 0.9. As mentioned, we have assumed an infinitesimally thin disk throughout this project. 
Accounting for the disk's finite thickness is necessary for understanding the behavior of de/dt for 
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e < 0.1. Second, as mentioned previously, the function de/dt may behave differently than we have 
calculated here for high eccentricities (e > 0.5). Our calculation will therefore benefit from an 
improved understanding of disk properties when a massive planet is present in a highly eccentric 
orbit. Third, additional disk effects, such as turbulence, may reduce the degree of corotation 
resonance saturation proposed by Ogilvie and Lubow to a level consistent with eccentricity growth 
on a timescale of a few million years, allowing planetary orbits to survive Type II migration. Finally, 
these results can be utilized in numerical studies of planet migration, such as MA05, where the 
effects of the disk can be combined with other solar system phenomena, such as secular resonances 
between planets. 

As this set of calculations illustrates, different disk properties, gap properties, and planet 
properties can lead to different types of behavior concerning the eccentricity evolution of giant 
planet orbits. In the vast ensemble of star and planet forming environments in our Galaxy (and 
others), we expect a wide variety of disk surface density profiles, gap widths, gap edge shapes, 
and other characteristics that determine eccentricity damping and excitation. In addition, for any 
given set of disk/planet properties, the eccentricity damping and/or excitation rates are complicated 
functions of eccentricity. As a result, although we have found a common denominator of eccentricity 
damping at e = 0.3 for a variety of gap shapes and widths, the question of how circumstellar disks 
affect the eccentricity evolution of their planets should generally be considered only in a statistical 
sense. 
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Coorbital Corotation 
Gap Material Saturation \i - m\ < 1 only all i, m 

Architecture Included? Included? 



wide and 
sharp-edged 

narrow and 
sharp-edged 

narrowest and 
sharp-edged 

Bate et al. gap 

Bate et al. gap 

narrow and 
sharp-edged 

narrowest and 
sharp-edged 

Bate et al. gap 

Bate et al. gap 



no 

yes 



no 

yes 



no 
no 

no 

no 
no 

yes 

yes 

yes 
yes 



-0.16 
-0.17 

-0.13 

-0.052 
-0.050 

0.063 

0.077 

0.032 
0.032 



-0.05 
-0.048 

-0.073 

-0.025 
-0.019 

0.029 

0.031 

0.017 
0.017 



Table 1: Here we present the timescale for growth or decay in eccentricity resulting from our calcula- 
tions. All numbers presented are in units of Myr, and negative values denote eccentricity damping 
while positive values denote eccentricity excitation. We present both the results of calculations 
including resonances with \£ — m\ < 1 and including all contributing resonances. We find that 
including resonances for which |^ — m| > 1 decreases the timescale by approximately a factor of 3, 
but does not change any cases from eccentricity damping to eccentricity excitation or vice versa. 
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Fig. 1. — The exact solution for (j)4^ ^{r\\_p>) (solid black curve) and an approximation for (1)4^ ^{r\\_p>) 
accurate to first order in eccentricity (dashed red line) as a function of eccentricity. This plot shows 
that, for this I, m, the linear approximation is useful only in the region e < 0.3. As the azimuthal 
wavenumber m increases, the cusp in the exact solution moves leftward while the oscillations in the 
right half of cl)^^{r\iR) increase in number and decrease in amplitude, thus shrinking the region of 
accuracy for the linear approximation. 
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Fig. 2. — The exact solution for the torque exerted by the inner Lindblad resonance Tf,4(riLR) as 
a function of e for i = 5, 6, 7, 8, 9, 10, 11, 12. As i increases, the first peak decreases and moves 
toward higher eccentricity. We find that by including resonances for which |£ — m| > 1, we are 
likely to obtain de/dt with increased magnitude, but are unlikely to significantly alter the shape of 
de/dt as a function of e. 
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Fig. 3. — The function de/dt versus e for a. 1 Mj planet in a 1 AU orbit. Here we have smoothed 
over the viscous scale length to demonstrate the shape of de/dt for a small smoothing length. The 
features of this function clearly result from the 4>^^ functions; each sharp peak and trough occurs 
at the location of a cusp in some 4>^^ function, which is the same location at which the radial 
derivative of 4>^^j^ becomes unbounded (e.g., see Fig. 1). 
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Fig. 4. — The function de/dt versus e for a 1 Mj planet in a 1 AU orbit. Here we have smoothed 
over the disk scale height to account for the finite thickness of the disk. The overall shape of this 
function is common for all surface density profiles that have gaps with sharp edges centered on the 
orbiting planet. The features of this function result from the 4>^^ functions; peaks and troughs 
for e < 0.5 occur at the locations of smoothed cusps in the 4>^^ functions, and result from the 
oscillations in low order c/)^^ functions for e > 0.5. The dotted red curve is for a narrower gap. 




Fig. 5. — Surface density profile for a gap as calculated in Bate et al. (2003). This profile is used 
to evaluate de/dt versus e in Fig. 6. 
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Fig. 7. — The function de/dt versus e for a 1 Mj planet in a 1 AU orbit for the gap shape of Fig. 
5 (from Bate et al. 2003). Here we have included coorbital corotation and Lindblad resonances as 
well as non-coorbital resonances. 
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Fig. 8. — Radial surface density profile for a disk in the presence of 1 Mj planets at (dashed blue) 
1 AU, (dotted red) 1.59 AU, and at both 1 AU and 1.59 AU (solid black). This configuration places 
the planets in a 2-1 mean motion resonance. The gaps produced by the two planets in the final 
case overlap to produce a single wide gap, in which each planet is offset from the gap center. 
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Fig. 9. — The resulting plot of de/dt as a function of e produced by the double gap shown in Fig. 
8. The dotted red line traces the eccentricity evolution of the outer planet, and the dashed blue 
line traces that of the inner planet. The function de/dtis changed substantially for the inner planet 
due to the lack of material at locations corresponding to the strongly contributing m = 2,i = 1 
Lindblad and corotation resonances. 
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Fig. 10. — The function de/dt versus e for a 1 Mj planet in a 1 AU orbit for the gap shape of Fig. 
5 (from Bate et al. 2003), including the effects of saturation of corotation resonances. Here we 
have included only non-coorbital corotation and Lindblad resonances. A marked difference exists 
between this graph and previous graphs. This is due to the reduction of non-coorbital corotation 
resonant torques via corotation saturation. 
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Fig. 11. — The function de/dt versus e for a 1 Mj planet in a 1 AU orbit for the gap shape of 
Fig. 5 (from Bate et al. 2003), including the effects of saturation of corotation resonances. Here 
we have included coorbital corotation and Lindblad resonances as well as non-coorbital resonances, 
although we then assume that coorbital corotation resonances are fully saturated. 
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Fig. 12. — Surface density as a function of radius for different gap shapes corresponding to values 
of eccentricity evolution timescales presented in Table 1 
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Fig. 13. — The function de/dt versus e for a 1 Mj planet in a 1 AU orbit for various disk configu- 
rations. Note that the range in eccentricity (horizontal axis) is restricted to the regime where the 
approximations of this paper are most applicable. The solid black curve corresponds to a gap in the 
disk with sharp edges defined by Tq = 70 K at ro = 7 AU, and the dotted red curve corresponds 
to a gap in the disk with sharp edges defined by Tq = 50 K at ro = 7 AU. The dashed green 
curve corresponds to the gap shape of Fig. 5 (from Bate et al. 2003) with no coorbital resonances 
contributing, and the long-dashed blue curve to the same disk shape with coorbital resonances 
contributing. In all four curves, corotation resonances are completely unsaturated. 



